# Packageslibrary(lubridate) # date supportlibrary(raster) # raster supportlibrary(rnaturalearth) # coastline polygonslibrary(sf) # simple feature supportlibrary(stars) # raster plotting with ggplotlibrary(gmRi) # styling for GMRIlibrary(scales) # axis labelslibrary(here) # project navigationlibrary(janitor) # data cleaninglibrary(patchwork) # multiplot arrangementlibrary(tidyverse) # data wrangling and plottinglibrary(xaringanExtra) # tab panelslibrary(ggthemes) # theme support# Support FunctionssuppressPackageStartupMessages(source(here("R/oisst_support_funs.R")))suppressPackageStartupMessages(source(here("R/temp_report_support.R")))path_fun <-os_fun_switch(mac_os ="mojave")# File Pathsmills_path <-path_fun("mills")oisst_path <-path_fun("res", "OISST/oisst_mainstays")# Set ggplot theme for figurestheme_set(theme_bw() +theme(# Titlesplot.title =element_text(hjust =0, face ="bold", size =14),plot.subtitle =element_text(size =9),plot.caption =element_text(size =7.2, margin =margin(t =20), color ="gray40"),legend.title =element_text(size =9),legend.text =element_text(size =7.5),# Axesaxis.line.y =element_line(),axis.ticks.y =element_line(), axis.line.x =element_line(),axis.ticks.x =element_line(), axis.text =element_text(size =11),# Facetsstrip.text =element_text(color ="white", face ="bold",size =11),strip.background =element_rect(color ="white", fill ="#00736D", size =1, linetype="solid") ) )# Polygons for mappingnew_england <-ne_states("united states of america", returnclass ="sf")canada <-ne_states("canada", returnclass ="sf")world <-ne_countries(returnclass ="sf")greenland <-ne_states(country ="greenland", returnclass ="sf")# Adding Logo to plotlogo_path <-paste0(system.file("stylesheets", package ="gmRi"), "/gmri_logo.png")lab_logo <- magick::image_read(logo_path)
Code
# Set parameters:params <-list(region ="apershing_gulf_of_maine",in_C =TRUE)# File paths for various extents based on params$regionregion_paths <-get_timeseries_paths(region_group ="gmri_sst_focal_areas", mac_os ="mojave")# clean up the nametidy_name <-str_replace_all(params$region, "_", " ") %>%str_to_title() %>%str_replace(pattern ="Of", replacement ="of")tidy_name <-str_remove(tidy_name, "Aak ")tidy_name <-str_remove(tidy_name, "Apershing ")tidy_name <-str_remove(tidy_name, "Cpr ")# Current Yearcurrent_yr <-year(Sys.Date())# Use panelsetuse_panelset()
Code
# Turn on the stylesheet:gmRi::use_gmri_style_rmd(css ="gmri_rmarkdown.css")
Code
# Turn on the gmri font for plotsshowtext::showtext_auto()
1: Global Temperature Maps
The 2022 global sea surface temperature anomalies have been loaded and displayed below to visualize how different areas of the ocean experience swings in temperature.
For perspective on where excess heat is distributed around the world, here first are maps of anomalies at a global perspective:
Code
# Access information to netcdf on boxnc_year <-"2022"anom_path <-str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", nc_year, ".nc")# Load 2022 as stackanoms_2022 <-stack(anom_path)
Annual Average
Code
# Make Annual Averageann_anom_ras <-calc(anoms_2022, fun = mean, na.rm = T)# Color scale titlesst_lab <-"Sea Surface Temperature Anomaly"# Color limit for palettestemp_limits <-c(-2, 2)temp_breaks <-c(temp_limits[1], temp_limits[1]/2, 0, temp_limits[2]/2, temp_limits[2])temp_labels <-str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")# Make st objectann_anom_st <-st_as_stars(rotate(ann_anom_ras)) # Build Mapggplot() +geom_stars(data = ann_anom_st) +geom_sf(data = world, fill ="gray30", color ="white", size = .25) +scale_fill_distiller(palette ="RdBu", na.value ="transparent", limit = temp_limits, oob = scales::squish,breaks = temp_breaks, labels = temp_labels) +guides("fill"=guide_colorbar(title = sst_lab, title.position ="right", title.hjust =0.5,barheight =unit(3, "inches"), frame.colour ="black", ticks.colour ="black")) +coord_sf(expand = F) +map_theme() +theme(legend.position ="right", legend.title =element_text(angle =90)) +labs(title =str_c("Global Temperature Anomalies - ", nc_year, " through: ", Sys.Date()))
Winter
Code
# Get previous year winterlast_nc_yr <-as.numeric(nc_year) -1last_yr_anom_path <-str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", last_nc_yr, ".nc")# Load previous year as stack to get the full winterlast_yr_anoms <-stack(last_yr_anom_path)# Join to 2022anoms_double <-stack(list(last_yr_anoms, anoms_2022))# Drop December 2022 so its not influencing the Winter of 2021-2022dec_2022 <-str_c("X", nc_year, ".12")not_dec21 <-which(str_sub(names(anoms_double), 1, 8) != dec_2022)anoms_nodec <- anoms_double[[not_dec21]]# Set up list of year and month combos to use map() for each seasonmonth_key <-list("Winter"=c(str_c(last_nc_yr, ".12"),str_c(nc_year, ".01"), str_c(nc_year, ".02")),"Spring"=str_c(nc_year, c("03", "04", "05"), sep ="."),"Summer"=str_c(nc_year, c("06", "07", "08"), sep ="."),"Fall"=str_c(nc_year, c("09", "10", "11"), sep ="."))
Code
# Get mean anoms across seasonsseason_stacks <-map(month_key, function(mon){# Get names from the stack - yyyy.mm stack_names <-names(anoms_nodec) stack_months <-str_sub(stack_names, 2,8)# layers with correct months: in_season <-which(stack_months %in% mon)# If there is no dates yet just multiply by zero so we can plot a blankif(length(in_season) ==0){ season_mean <- anoms_nodec[[1]] *0 } else {# Get mean across those months season_mean <-calc(anoms_nodec[[in_season]], mean, na.rm = T) }# season_mean <- st_as_stars(rotate(season_mean))return(season_mean)})
The regional views of 2022’s sea surface temperature anomalies have been loaded and displayed below to visualize how localized differences changed throughout the year.
Code
# Load the bounding box for Andy's GOM to show they alignpoly_path <- region_paths[[params$region]][["shape_path"]]region_extent <-st_read(poly_path, quiet =TRUE)# Pull extents for the region to set crop extentcrop_x <-st_bbox(region_extent)[c(1,3)]crop_y <-st_bbox(region_extent)[c(2,4)]# # Expand the area out to see the larger patternscrop_x <- crop_x +c(-2.5, 3.5)crop_y <- crop_y +c(-1.5, 1)# Make a new extent for croppingregion_extent_expanded <-st_sfc(st_polygon(list(rbind(c(crop_x[[1]], crop_y[[1]]), c(crop_x[[1]], crop_y[[2]]), c(crop_x[[2]], crop_y[[2]]), c(crop_x[[2]], crop_y[[1]]), c(crop_x[[1]], crop_y[[1]])))))# convert to sfregion_extent_expanded <-st_as_sf(region_extent_expanded)
Looking specifically at the last heatwave event, we can step through how the event progressed over time, and developing pockets or warmer/colder water masses.
Code
# Identify the last heatwave event that happenedlast_event <-max(region_hw$mhw_event_no, na.rm = T)# Pull the dates of the most recent heatwavelast_event_dates <- region_hw %>%filter(mhw_event_no == last_event) %>%pull(time)# Buffer the dates, start 7 days beforeevent_start <- (min(last_event_dates) -7)event_stop <-max(last_event_dates)date_seq <-seq.Date(from = event_start,to = event_stop,by =1)# Load the heatwave datesdata_window <-data.frame(time =c(min(date_seq) , max(date_seq) ),lon = crop_x,lat = crop_y)# Pull data off boxhw_stack <-oisst_window_load(data_window = data_window, anomalies = T, mac_os ="mojave")#drop any empty years that bug inhw_stack <- hw_stack[map(hw_stack, class) !="character"]##### Format the layers and loop through the maps ##### Grab only current year, format datesthis_yr <-stack(hw_stack)day_count <-length(names(this_yr))day_labs <-str_replace_all(names(this_yr), "[.]","-")day_labs <-str_replace_all(day_labs, "X", "")day_count <-c(1:day_count) %>%setNames(day_labs)# Progress through daily timeline to indicate heatwave status and severityhw_timeline <- region_hw %>%filter(time %in%as.Date(day_labs))
Code
#### Plot Settings:# Set palette limits to center it on 0 with scale_fill_distillerlimit <-c(max(values(this_yr), na.rm = T) *-1, max(values(this_yr), na.rm = T) )# Plot Heatwave 1 day at a time as a GIFday_plots <-imap(day_count, function(date_index, date_label) {# grab dates heatwaves_st <-st_as_stars(this_yr[[date_index]])#### 1. Map the Anomalies in Space day_plot <-ggplot() +geom_stars(data = heatwaves_st) +geom_sf(data = new_england, fill ="gray30", color ="white", size = .25) +geom_sf(data = canada, fill ="gray30", color ="white", size = .25) +geom_sf(data = greenland, fill ="gray30", color ="white", size = .25) +geom_sf(data = region_extent, color =gmri_cols("gmri blue"), linetype =2, size =1,fill ="transparent") +scale_fill_distiller(palette ="RdYlBu", na.value ="transparent", limit = limit) +map_theme(legend.position ="bottom") +coord_sf(xlim = crop_x, ylim = crop_y, expand = T) +guides("fill"=guide_colorbar(title ="Sea Surface Temperature Anomaly \u00b0C",title.position ="top",title.hjust =0.5,barwidth =unit(3, "in"),frame.colour ="black",ticks.colour ="black")) # Set colors by name color_vals <-c("Sea Surface Temperature"="royalblue","Heatwave Event"="darkred","MHW Threshold"="coral3","Daily Climatology"="gray30")#### 2. Plot the day and the overall anomaly to track dates date_timeline <-ggplot(data = hw_timeline, aes(x = time)) +geom_line(aes(y = sst, color ="Sea Surface Temperature")) +geom_line(aes(y = hwe, color ="Heatwave Event")) +geom_line(aes(y = cse, color ="Cold Spell Event")) +geom_line(aes(y = mhw_thresh, color ="MHW Threshold"), lty =3, size = .5) +geom_line(aes(y = mcs_thresh, color ="MCS Threshold"), lty =3, size = .5) +geom_line(aes(y = seas, color ="Daily Climatology"), lty =2, size = .5) +scale_color_manual(values = color_vals) +# Animated Point / linegeom_point(data =filter(hw_timeline, time ==as.Date(date_label)),aes(time, sst, shape =factor(mhw_event)), color =gmri_cols("gmri blue"), size =3, show.legend =FALSE) +geom_vline(data =filter(hw_timeline, time ==as.Date(date_label)),aes(xintercept = time), color ="gray50", size =0.5,linetype =3,alpha =0.8) +labs(x ="", y ="",color ="",subtitle ="Regional Temperature \u00b0C",shape ="Heatwave Event") +theme(legend.position ="bottom")#### 3. Assemble plot(s) p_layout <-c(area(t =1, l =1, b =2, r =8),area(t =3, l =1, b =8, r =8))# plot_agg <- (date_timeline / day_plot) + plot_layout(heights = c(1, 3)) plot_agg <- date_timeline + day_plot +plot_layout(design = p_layout)return(plot_agg )})walk(day_plots, print)
3: Ranking the Rate of Warming
If we look at the rates of change from 1982-2021 for each grid cell, rather than the observed temperature, it is possible to rank how hot each location on earth is warming relative to all the others.
Once we have the rankings, we can then take the average ranking within the Gulf of Maine we can obtain the average warming rank for the area compared to the rest of the globe.
# Get the rank information that go with the original extent used by the timelinesranks_masked <-mask_nc(ranks_stack_all, region_extent)rates_masked <-mask_nc(rates_stack_all, region_extent)region_ranks <-get_masked_vals(ranks_masked, rates_masked)# Prep it for text input.avg_rank <- region_ranks$`Mean Rank`*100avg_rate <- region_ranks$`Mean Rate`low_rank <- region_ranks$`Min Rank`*100low_rate <- region_ranks$`Min Rate`top_rank <- region_ranks$`Max Rank`*100top_rank <-ifelse(top_rank ==100, "greater than or equal to 99.5", top_rank)top_rate <- region_ranks$`Max Rate`
Based on data from 1982-2021, the warming rates of Gulf of Maine have been some of the highest in the world. The area as a whole has been increasing at a rate of 0.044\(^{\circ}C/year\) which is faster than 95.9% of the world’s oceans.
Over that same period locations within the Gulf of Maine have been warming at rates as low as 0.017\(^{\circ}C/year\) and as rapidly as 0.094\(^{\circ}C/year\), corresponding to ranks as low as 60.3% and as high as greater than or equal to 99.5%.
Mapped below are the corresponding warming rates and their global rankings.
Code
# For the full map we mask again, but zoom out a little# Mask again using the expanded mask so we can zoom outranks_masked <-mask_nc(ranks_stack_all, region_extent_expanded)rates_masked <-mask_nc(rates_stack_all, region_extent_expanded)# Make stars objectrank_stars <-st_as_stars(ranks_masked)rates_stars <-st_as_stars(rates_masked)
